simulate_RAR <- function(N, K, mu_1, mu_2, sigma, I_theta_fix,delta) {
  # Ratio, estimates of mu, theta
  muhat_1 <- muhat_2 <- ratio <- numeric(K)
  theta_hat <- rep(NA, K)
  # data from trial
  x_1 <- x_2 <- numeric(0)
  # Information
  I1 <- I2 <- numeric(K)
  # Number of patients per arm
  n1 <- n2 <- n1.new <- n2.new <- w <- numeric(K)
  S1 <- N1 <- S2 <- N2 <- 0
  for (k in 1:K) {
    if (k == 1) {
      # Initial equal allocation
      ratio[k] <- 1
      # Enroll 20 patients when N = 100, K=5
      n1.new[k] <- n2.new[k] <- N / K
      n1[k] <- n2[k] <- N / K
    } else {
      # Define ratio based on formula with previous estimate of theta
      ratio[k] <-  a ** (theta_hat[k - 1]/ (2 * delta))
      # Patients in N1
      n1[k] <- (k / 5) * (sigma ^ 2) * I_theta_fix * (1 + ratio[k])
      # Patient in N2
      n2[k] <- (k / 5) * (sigma ^ 2) * I_theta_fix * (1 + (1 / ratio[k]))
      # Difference in patients from previous
      n1.new[k] <- n1[k] - n1[k - 1]
      n2.new[k] <- n2[k] - n2[k - 1]
      if(n1.new[k]<=0){
        n1[k]= n1[k - 1]+0.00001
        n1.new[k] = n1[k] - n1[k - 1]
        #n2[k] = (((k/5)*I_theta_fix)^(-1)-1/(n2/sigma^2))^(-1)
        n2[k] = (k*n1[k]*sigma^2*I_theta_fix)/(5*n1[k]-k*sigma^2*I_theta_fix)
        n2.new[k] = n2[k]-n2[k-1]
      }
      if(n2.new[k]<=0){
        n2[k]= n2[k - 1]+0.00001
        n2.new[k] = n2[k] - n2[k - 1]
        #n1[k] = (((k/5)*I_theta_fix)^(-1)-1/(n2[k]/sigma^2))^(-1)
        n1[k] = (k*n2[k]*sigma^2*I_theta_fix)/(5*n2[k]-k*sigma^2*I_theta_fix)
        n1.new[k] = n1[k]-n1[k-1]
      }
    }
    #Sample once from normal distribution
    x_1[k] <- rnorm(1, mu_1, sqrt(sigma^2/n1.new[k]))
    x_2[k] <- rnorm(1, mu_2, sqrt(sigma^2/n2.new[k]))
    
    muhat_1[k] <- sum(n1.new[1:k] * x_1[1:k]) / sum(n1.new)
    muhat_2[k] <- sum(n2.new[1:k] * x_2[1:k]) / sum(n2.new)
    theta_hat[k] <- muhat_1[k] - muhat_2[k]
  }
  return(list(n1 = n1, n2 = n2, theta_hat = theta_hat))
}

posterior_mean <- function(prior_mean, prior_tau, data, tau, n,k) {
  result <- tryCatch({
    # Check if the lengths of n and data are compatible
    if (length(n) != length(data)) {
      stop("Length of n and data must be the same. Found n of length ", length(n), " and data of length ", length(data))
    }
    # Calculate post mean
    posterior_mean_result <- ((prior_tau) * prior_mean + sum(n * (tau) * (data))) / (sum(n* tau) + prior_tau)
    return(posterior_mean_result)
    
  }, error = function(e) {
    # Catch the error and return a message
    cat("Error in posterior_mean (K=):",k, e$message, "\n")
    cat("n: ",paste(n),"\n")
    cat("data: ",paste(data),"\n")
    return(NA)  # Return NA or any default value in case of error
  })
  
  return(result)
}

x_1[k] <- rnorm(1, mean = mu_1, sd = sqrt(sigma ^ 2 / n1.new[k]))
x_2[k] <- rnorm(1, mean = mu_2, sd =  sqrt(sigma ^ 2 / n2.new[k]))

# Posterior variances
tau_n1[k] <- posterior_tau(prior_tau = tau_01,
                           tau = tau,
                           n = n1.new[1:k])

tau_n2[k] <- posterior_tau(prior_tau = tau_02,
                           tau = tau,
                           n = n2.new[1:k])
# Posterior means
mu_n1[k] <- posterior_mean(
  prior_mean = mu_01,
  n = n1.new[1:k],
  tau = tau,
  prior_tau = tau_01,
  data = x_1[1:k],
  k = k
)

mu_n2[k] <- posterior_mean(
  prior_mean = mu_02,
  n = n2.new[1:k],
  tau = tau,
  prior_tau = tau_02,
  data = x_2[1:k],
  k = k
)

theta_hat[k] <- mu_n1[k] - mu_n2[k]